{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 使用 Pyomo 求解背包问题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "【动态规划问题5】 某货运公司负责四种货物的运输任务，装载工具是一种最大承载能力为15吨的卡车，每种货物的质量及价值如表所列。\n",
    "\n",
    "货物编号 | 1 | 2 | 3 | 4 | 5\n",
    "-------- | --- | --- | --- | --- | ---\n",
    "单位质量 / 吨 | 3 | 4 | 2 | 5 | 4\n",
    "单位价值 / 百元 | 5 | 3 | 4 | 5 | 6\n",
    "\n",
    "请问：\n",
    "\n",
    "1. 试问使总价值最大的装载方案?\n",
    "2. 如果货物 1 和货物 2 至少装一件，试问使总价值最大的装载方案?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 导入相关的库\n",
    "\n",
    "需要 Pyomo 工具包。\n",
    "\n",
    "Pyomo 是一个基于 Python 的开源软件包，它支持多种优化功能，用于制定和分析优化模型。可以通过下面的命令安装：\n",
    "\n",
    "```bash\n",
    "pip install pyomo\n",
    "```\n",
    "\n",
    "此处以 [SCIP 求解器](https://scipopt.org/index.php#download) 为例。实际上对于 Pyomo 而言，各求解器可通用。 **求解器需要单独下载、安装并将安装路径添加到环境变量。**\n",
    "\n",
    "通过 `pyomo.environ` 导入 Pyomo 建模描述语言环境。\n",
    "\n",
    "`SolverFactory` 可以允许我们选择和使用特定的求解器。无需生成 `.lp` 文件，而是可以直接在 Python 编程环境中调用求解器。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyomo.environ import *\n",
    "from pyomo.opt import SolverFactory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第一问\n",
    "\n",
    "设置五个整数规划变量 $x_i \\in \\mathbb{Z^+}, i \\in \\{1, 2, 3, 4, 5\\}$，分别表示每一种货物装载的件数。$v_i$ 表示每一种货物的价值，$w_i$ 表示每一种货物的重量。则原问题可以用如下的规划模型描述：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max{ f^{*} } =& \\sum_{i=1}^{5}{(v_i \\times x_i)} \\\\\n",
    "\\mathrm{s.t. \\ \\ \\ \\ }\n",
    "& \\sum_{i=1}^{5}{(w_i \\times x_i)} \\le C \\\\\n",
    "& x_i \\in \\mathbb{Z^{+}}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "编写 Pyomo 代码求解："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先，针对问题创建 Pyomo 规划问题模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ConcreteModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "设置规划问题的各个参数。其中：\n",
    "\n",
    "- `N` 为模型中规划变量的总个数\n",
    "- `I` 为一个 `range(1, 5)` 对象，用于在 `pyomo.Var()` 方法中一次性创建 5 个规划变量\n",
    "- `C` 为总容量\n",
    "- `W` 和 `v` 即原始问题中的货物的单位质量和单位价值。\n",
    "\n",
    "这些参数理论上也可以通过 `Param()` 等方法实例化为模型 `model` 的参数（属性）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5\n",
    "I = range(1, N+1)\n",
    "\n",
    "C = 15\n",
    "\n",
    "W = [3, 4, 2, 5, 4]\n",
    "v = [5, 3, 4, 5, 6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为 `model` 增添一组决策变量 `x`，使用 `I` 创建。可以通过 `model.x[i]` 这样的形式来索引，但是查看 `x` 的值则需要使用到 `value()` 函数，即：`value(model.x)`。\n",
    "\n",
    "同时，设置变量类型为 `NonNegativeIntegers`，表示规划变量是一组非负整数变量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.x = Var(\n",
    "    I, bounds=(0, 15), \n",
    "    domain=NonNegativeIntegers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为模型添加目标函数。`expr` 参数表示“表达式”。此处可以通过在 `sum()` 中使用列表推导式来表示 $\\sum$ 求和的含义。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.obj = Objective(\n",
    "    expr=(\n",
    "        sum(model.x[i] * v[i-1] for i in I)\n",
    "        ), sense=maximize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为模型添加约束条件。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.con = Constraint(\n",
    "    expr = (\n",
    "        sum(model.x[i] * W[i-1] for i in I) <= C\n",
    "        ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "调用 SCIP 求解器求解。返回的 `results` 可以查看到模型求解的输出结果。\n",
    "\n",
    "可以看到，模型的返回信息是 `optimal solution found`，其中，找到的最优解的值 $f^*$ 为 29。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Problem: \n",
      "- Lower bound: -inf\n",
      "  Upper bound: 29.0\n",
      "  Number of objectives: 1\n",
      "  Number of constraints: 0\n",
      "  Number of variables: 5\n",
      "  Sense: unknown\n",
      "Solver: \n",
      "- Status: ok\n",
      "  Message: optimal solution found\n",
      "  Termination condition: optimal\n",
      "  Id: 0\n",
      "  Error rc: 0\n",
      "  Time: 0.0\n",
      "  Gap: 0.0\n",
      "  Primal bound: 29.0\n",
      "  Dual bound: 29.0\n",
      "Solution: \n",
      "- number of solutions: 0\n",
      "  number of solutions displayed: 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "opt = SolverFactory('scip') # 调用 SCIP\n",
    "results = opt.solve(model) # 求解模型\n",
    "\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面的代码可以输出模型的参数信息，查看模型的变量、目标和约束的声明情况，以及对应的当前值。\n",
    "\n",
    "由于已经求解了模型，模型求解的值也会保存在 `model` 当中，所以这个时候我们输出模型的信息，可以看见各个规划变量的值。从输出结果来看，求解器为我们找到的规划方案是装载一件商品 1 和 6 件商品 3。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Set Declarations\n",
      "    x_index : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    5 : {1, 2, 3, 4, 5}\n",
      "\n",
      "1 Var Declarations\n",
      "    x : Size=5, Index=x_index\n",
      "        Key : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "          1 :     0 :   1.0 :    15 : False : False : NonNegativeIntegers\n",
      "          2 :     0 :  -0.0 :    15 : False : False : NonNegativeIntegers\n",
      "          3 :     0 :   6.0 :    15 : False : False : NonNegativeIntegers\n",
      "          4 :     0 :  -0.0 :    15 : False : False : NonNegativeIntegers\n",
      "          5 :     0 :   0.0 :    15 : False : False : NonNegativeIntegers\n",
      "\n",
      "1 Objective Declarations\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : maximize : 5*x[1] + 3*x[2] + 4*x[3] + 5*x[4] + 6*x[5]\n",
      "\n",
      "1 Constraint Declarations\n",
      "    con : Size=1, Index=None, Active=True\n",
      "        Key  : Lower : Body                                       : Upper : Active\n",
      "        None :  -Inf : 3*x[1] + 4*x[2] + 2*x[3] + 5*x[4] + 4*x[5] :  15.0 :   True\n",
      "\n",
      "4 Declarations: x_index x obj con\n"
     ]
    }
   ],
   "source": [
    "model.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第二问\n",
    "\n",
    "第二问要求货物 1 和货物 2 至少装一件，因此在原始的规划问题的基础上添加两条约束，限制第一和第二种货物的数量大于等于 1。\n",
    "\n",
    "针对第二问建立的规划模型如下：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max{ f^{*} } =& \\sum_{i=1}^{5}{(v_i \\times x_i)} \\\\\n",
    "\\mathrm{s.t. \\ \\ \\ \\ }\n",
    "& \\sum_{i=1}^{5}{(w_i \\times x_i)} \\le C \\\\\n",
    "& x_1 \\ge 1 \\\\\n",
    "& x_2 \\ge 1 \\\\\n",
    "& x_i \\in \\mathbb{Z^{+}}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "初始化规划模型、变量等的过程与第一问相同。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ConcreteModel()\n",
    "model.x = Var(I, bounds=(0, 15), \n",
    "    domain=NonNegativeIntegers)\n",
    "model.obj = Objective(\n",
    "    expr=(sum(model.x[i] * v[i-1] for i in I)\n",
    "        ), sense=maximize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在模型 1 的基础上，增加两个约束。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.con = Constraint(\n",
    "    expr=(\n",
    "        sum(model.x[i] * W[i-1] for i in I) <= C\n",
    "    ))\n",
    "model.con_x1 = Constraint(\n",
    "    expr=(model.x[1] >= 1 ))\n",
    "model.con_x2 = Constraint(\n",
    "    expr=(model.x[2] >= 1 ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用求解器进行求解的结果如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Problem: \n",
      "- Lower bound: -inf\n",
      "  Upper bound: 24.0\n",
      "  Number of objectives: 1\n",
      "  Number of constraints: 0\n",
      "  Number of variables: 5\n",
      "  Sense: unknown\n",
      "Solver: \n",
      "- Status: ok\n",
      "  Message: optimal solution found\n",
      "  Termination condition: optimal\n",
      "  Id: 0\n",
      "  Error rc: 0\n",
      "  Time: 0.0\n",
      "  Gap: 0.0\n",
      "  Primal bound: 24.0\n",
      "  Dual bound: 24.0\n",
      "Solution: \n",
      "- number of solutions: 0\n",
      "  number of solutions displayed: 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "opt = SolverFactory('scip') # 调用 SCIP\n",
    "results = opt.solve(model) # 求解模型\n",
    "\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从求解后的模型结果来看，这里的决策是运输了第一、第二种商品各一件；第三种商品四件，最大价值 $f^*$ 为 24。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 Set Declarations\n",
      "    x_index : Size=1, Index=None, Ordered=Insertion\n",
      "        Key  : Dimen : Domain : Size : Members\n",
      "        None :     1 :    Any :    5 : {1, 2, 3, 4, 5}\n",
      "\n",
      "1 Var Declarations\n",
      "    x : Size=5, Index=x_index\n",
      "        Key : Lower : Value : Upper : Fixed : Stale : Domain\n",
      "          1 :     0 :   1.0 :    15 : False : False : NonNegativeIntegers\n",
      "          2 :     0 :   1.0 :    15 : False : False : NonNegativeIntegers\n",
      "          3 :     0 :   4.0 :    15 : False : False : NonNegativeIntegers\n",
      "          4 :     0 :  -0.0 :    15 : False : False : NonNegativeIntegers\n",
      "          5 :     0 :   0.0 :    15 : False : False : NonNegativeIntegers\n",
      "\n",
      "1 Objective Declarations\n",
      "    obj : Size=1, Index=None, Active=True\n",
      "        Key  : Active : Sense    : Expression\n",
      "        None :   True : maximize : 5*x[1] + 3*x[2] + 4*x[3] + 5*x[4] + 6*x[5]\n",
      "\n",
      "3 Constraint Declarations\n",
      "    con : Size=1, Index=None, Active=True\n",
      "        Key  : Lower : Body                                       : Upper : Active\n",
      "        None :  -Inf : 3*x[1] + 4*x[2] + 2*x[3] + 5*x[4] + 4*x[5] :  15.0 :   True\n",
      "    con_x1 : Size=1, Index=None, Active=True\n",
      "        Key  : Lower : Body : Upper : Active\n",
      "        None :   1.0 : x[1] :  +Inf :   True\n",
      "    con_x2 : Size=1, Index=None, Active=True\n",
      "        Key  : Lower : Body : Upper : Active\n",
      "        None :   1.0 : x[2] :  +Inf :   True\n",
      "\n",
      "6 Declarations: x_index x obj con con_x1 con_x2\n"
     ]
    }
   ],
   "source": [
    "model.pprint()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
